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1. Introduction 

Among the various analytical and numerical approaches to strongly correlated systems, 
numerical exact diagonalization takes a unique position. It is not burdened by any 
assumptions or approximations and thus provides unbiased benchmarks for other 
analytical and/or numerical approaches [1]. It is also appealing in its conceptually 
simple and straightforward nature. The basic idea, to set up the Hamiltonian matrix in 
some basis and thus reduce a physical problem to a purely mathematical one, is readily 
accessible to a senior undergraduate student. 

However, possibly due to some technical subtleties, exact diagonalization is not 
accounted for in detail in existing textbooks on computational physics. It is the 
aim of this paper to illustrate these tricks and promote teaching and using of exact 
diagonalization. To make the discussion concrete, we take the Bose-Hubbard model as 
an example. This model is chosen because of its relevance to the currently active field of 
ultracold atom physics [2] . It has been realized with ultracold atoms in an optical lattice 
and the celebrated Superfiuidity-Mott insulator (SF-MI) transition has been observed 
experimentally [3l H] . We will use exact diagonalization to get a glimpse of this quantum 
phase transition. 
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One common misconception, according to the experience of the authors, is that in 
doing numerical exact diagonahzation, one solves all the eigenvalues and eigenvectors 
of the Hamiltonian by some algorithm. This ideal case is actually neither possible nor 
necessary in many cases as long as the dimension of the Hilbert space D gets large. It 
is impossible since to reduce a Hermitian matrix H in the form UAW , with A being a 
diagonal real matrix and U a unitary matrix, it would take time on the order of 0{D^) 
and memory space on the order of 0{D'^). With a moderate value D = 100 000, the 
memory needed is over 10 GB, far beyond that of a typical desktop computer, needless 
to say the time cost. It is also unnecessary since physically, in many cases, the most 
relevant eigenstates are the ground state and low lying excited states. High excited 
states, due to the Boltzmann factor, contribute little to the thermodynamics of the 
system in low temperatures. 

In view of the considerations mentioned above, one can fully appreciate the value 
of the Lanczos algorithm |5]. This algorithm belongs to the iterative category for 
solving eigenvalue problems. As the iteration goes on, the estimated eigenvalues and 
eigenvectors converge quickly. Especially, the extremal eigenvalues and eigenvectors 
converge first. Usually, with an iteration time m D, the ground state and several low 
excited states converge to machine precision. In a certain sense, the Lanczos algorithm 
is a tailor-made algorithm for solving the ground state and/or low lying excited states 
of a Hamiltonian. It provides exactly what we need for us, no more no less. 

As far as we know, all exact diagonalizations are based on the Lanczos algorithm 
and its variants. Since this algorithm has become a standard topic in textbooks on 
numerical matrix theory [6] and since there are many monographes [7] devoted to this 
algorithm and also several very readable introductions [HI E], here in this article, we 
would not go into the details of this algorithm. We will just invoke some packages based 
on Lanczos algorithm and use it as the final stroke. 

On the contrary, our emphasis is placed on some other techniques which we believe 
are involved in all kinds of exact diagonalizations in a wide variety of contexts. The 
road map we will take is perhaps the most natural one — first enumerate all the basis 
vectors, then set up the Hamiltonian matrix in this basis, and finally invoke the Lanczos 
algorithm to solve the desired eigenstates. In each step, we will explain the tricks in 
detail. We would like to mention that in practice, usually the Hamiltonian matrix is not 
explicitly set up beforehand, instead the action of the Hamiltonian on a wave vector is 
done "on the fly" [HI E] . This is consistent with the philosophy of iterative method — the 
matrix-vector multiplication is all what we need and its internal workings are of no 
concern. Though this is of use for saving the memory, we would not introduce it here 
for our pedagogical purposes. 
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2. The model and its symmetries 

We begin by describing tlie one dimensional Bose-Hubbard model and its symmetries. 
The Hamiltonian is 

U 

H = - J ^(ajflj + ajflj) + — ^ni{ni - 1) (1) 

(ij) «=1 

where a] (a,) creates (annihilates) a particle on the site i and hi = a\ai counts the 
particle number on that site. The first term proportional to J is the kinetic part of the 
Hamiltonian (Hkin) and describes particle hopping between adjacent sites (in the sum, 
(ij) = {ji)). The second term is the interaction part (Hint) and is due to the particle- 
particle interaction, the strength of which is characterized by the parameter U. The 
Bose-Hubbard model has been realized with ultracold boson atoms in an optical lattice 
[1]. Moreover, in this system, the parameters J and U can be conveniently adjusted by 
various means, e. g. the Feshbach resonance or just changing the intensity of the laser 
beams. 

The Hamiltonian H possesses several symmetries. The first one is the U{1) 
symmetry, which is associated with the conservation of the total atom number N = 
YlfLi'f^i- The Hamiltonian is invariant under the transform (aj,ai) — )■ (a|e*^, ajC"*^) = 
e*^^(aj, ai)e~*^^ for V ^ G M. The second one is the translation symmetry. The 
Hamiltonian is invariant under the transform (a|,aj) — >■ (a|^^,aj+i), if the periodic 
boundary condition is imposed. This symmetry is associated with the conservation of 
the total quasi-momentum of the system 

^^°^2vr), (2) 

q=0 ^ ^ 

where the operator 

1 

fet = ^ye^0-2-^A^)a| (3) 

creates a particle in the Bloch state with quasi-momentum 27rg/M. Actually, the 
transform above is done as e~*^(a|, ai)e*^ = (aj_,_]^, aj+i). In terms of (6j,6q), the 
Hamiltonian H is rewritten as 

M-1 jj M-1 M-1 

H = -2jY,cos{27,q/M)h\h, + — hlhlKKA.+,2,q,+q,, (4) 

q=Q gi,g2=0 (J3,ij4=0 

where the Dirac function is defined as 

_ 1 if gi + g2 = gs + ^4 (mod M) 
(J otherwise 

It is then clear that K is conserved. The third symmetry is the reflection symmetry. 
The Hamiltonian is also invariant under the transform (a|,aj) — )■ (aj^,/_j, ctM-j) [10] - 
or in terms of (6j,6g), {h\,hq) — > (6l.g,6_g). Combination of the translation and 
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reflection symmetries indicates that the Bose-Hubbard model has the Dm symmetry, 
the symmetry of a equilateral polygon with M vertices. This is plausible if we envisage 
that the M sites are placed equidistantly on a circle. 

Therefore, the Bose-Hubbard model is of U{1) ® Dm symmetry. It is desirable 
to decompose the total Hilbert space into subspaces according to the irreducible 
representations of this group. The Hamiltonian cannot couple two states belonging 
to two different irreducible representations and thus is block (partially) diagonalized 
P, [H]. Analytically, it can be proven that the ground state \G) of H belongs to the 
identity representation of Dm [12j. Therefore, as far as the ground state is concerned, 
we only need to seek it in a subspace where all the basis vectors are of a definite atom 
number [thus belong to a definite representation of U{1)] and are invariant under all 
rotations and reflections. 

We first restrict to the space "H with total atom number being N. The dimension 
of this space is found to be 

m{M-i)\ ' ^ ^ 

which grows explosively with the system size. For fixed filling factor N/M = 1, 
D = 24 310 for M = 9, and it grows to D = 352 716 for M = 11, and further to 
D = 5 200 300 for M = 13. We may divide this space into M smaller subspaces 
according to the eigenvalues of K. The ground state, being translationally invariant, falls 
in the subspace Ho with K = 0, whose dimension Do is approximately D/M. Actually, 
Do = 2 704, 32 066, and 400 024 in the case of = M = 9, 11, and 13, respectively (see 
reference [Mj)- The subspace Ho can be further divided into two subspaces according to 
the two representations of the reflection group {/, a}. The ground state, being invariant 
under reflection, belongs to the subspace T-Lq, where the superscript means all the basis 
vectors yield a plus sign under reflection. The dimension Dq of is nearly half of 
Dq. Actually, Dq = 1 387, 16 159, and 200 474 respectively in the three cases above. 
The reduction of the dimension from D to Do and again to D^ promises a reduction of 
computation, especially, a reduction of memory needed. 

Indeed, when memory is limited, it is necessary to work within the subspace Hq 
(or even T-Lq) and with the Hamiltonian in equation (j4]) (so we do in the A^ = M = 13 
case below). However, to simplify the discussion and focus on essential techniques, we 
will still work within the space "H and with the Hamiltonian in equation ([T]). In fact, 
working in the subspaces Ho or Hq requires a bit more effort in coding and will be left 
as exercises. 

3. Basis vectors generation 

A natural basis is the occupation number basis {|ni,n2, . . . jUm)} which are defined as 
hi\ni, ria, . . . , Um) = ni\ni, ns, . . . , Um) (7) 
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with Tii > 0. In the subspace with a fixed total particle number A^, we have the 
constraint: X^i^i — ^- need to enumerate all the basis vectors satisfying this 
constraint. One naive idea is to write down a piece of code with a M-fold loop: 

for n_l=0:N 

for n_2=0:N-n_l 

for n_3=0:N-n_l-n_2 

end 
end 
end 

This approach, though workable, has two apparent drawbacks. First, the number of 
loops depends on the number of sites and hence the code is infiexible. Second, when 
coding with tools such as MATLAB which is inefficient in dealing with loops, the 
efficiency would be low. 

Here we prescribe one way to bypass these difficulties. To this end, we first note 
that it is possible to rank all the basis vectors |ni,n2, . . . jUm) in lexicographic order 
[13]. For two different basis vectors |?t,i, n2, . . . , um) and \ni,n2, . . . ,nM) , there must 
exist a certain index 1 < k < M — 1 such that Ui = fii for 1 < i < k — 1 while ^ fik. 
We say n2, . . . , Um) is superior (inferior) to |ni, ^2, . . . , Um) if > Uk {n^ < n^). It 
can be shown that this defines a total order among the basis vectors. In particular, it is 
clear that |A^, 0, . . . , 0) is superior to all other basis vectors while |0, 0, . . . , A^) is inferior 
to all other basis vectors. 

Having furnished the set of basis vectors with an order structure, we can now 
generate all the basis vectors one by one by descending from the highest one |A^, 0, . . . , 0). 
Given a basis vector \ni,n2, . . . ,nM) with um < N, we proceed to the next basis vector 
inferior to the current one according to the following rule [H]: 

Suppose nfc 7^ while = for all /c + 1 < i < M — 1, then the next basis vector 
is \ni, 712, . . . , Um) with 

• fii = Hi for 1 < 2 < — 1; 

• nfc = nfc - 1; 

• fik+i = N — Yli=i ^^'i = for i>k + 2. 

This procedure will end with the lowest basis vector |0, 0, . . . , A^). Obviously, this 
algorithm yields a code involving only a single loop and with all the difficulties associated 
with the naive one avoided. Numerically, we store the basis vectors in a D x M array 
A, with the v-th generated basis vector filled in the f-th row of the array. We will refer 
to the basis vector in v-th. row as |f ), so \v) = \Ayi, A^2-, ■ ■ ■ , A^jm)- 

As an example, we enumerate in table [1] all the basis vectors generated with the 
foregoing algorithm in the case of A^ = M = 3. As for the efficiency of the algorithm, 
we mention that in the case of A^ = M = 13, it takes about 38 seconds to generate the 
D = 5 200 300 basis vectors with our MATLAB code in our desktop computer |15j . 
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Table 1. Configurations of the basis vectors 1 711,712,713) witli atom number iV = 3 
and site number AI — 3. Tliey are generated recursively according to the algorithm 
described in section [3] 
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4. Setting up the Hamiltonian matrix 

With all the basis vectors prepared, we are now in the position to set up the Hamiltonian 
matrix with respect to this basis. That is, we are to determine the D x D matrix H 
corresponding to the Bose-Hubbard Hamiltonian H with 

= {u\H\v). (8) 

Here by determining the matrix H, we do not mean to save it in the full matrix form 
in the computer (that will cost memory on the order of D^), but to figure out all its 
non-zero elements and their positions, i.e., their row and column numbers. Actually, 
as we will see below, the matrix H is extremely sparse with at most 2M + 1 non-zero 
elements per column. Therefore, it is appropriate to store H in a certain sparse form, 
which will require memory only on the order of D. In MATLAB, a sparse matrix is 
stored in the coordinate format. 

To proceed, we treat the interaction part Hint and kinetic part Hkin of the 
Hamiltonian separately. The corresponding matrices are denoted as H_int and H_kin, 
respectively. We note that H_int and H_kin are the diagonal and off-diagonal parts of 
H respectively. We also note that this separation is necessary when we want to change 
the ratio U/J to study the SF-MI transition. The matrix H_int can be easily done, 
therefore, we will concentrate on H_kin. 

A general and straightforward but naive method to set up H_kin is to let u and 
V run over all the integers from 1 to D, respectively, and examine the corresponding 
matrix elements one by one. This procedure entails computation scale proportional to 
D^, and is very inefficient since most checks yield null results. A clever way out is to ask 
the question, in each column, which elements are non-zero? Physically, it is equivalent 
to ask, given an arbitrary basis vector \v), if we act Hkin on it, which (generally not 
merely one) basis vectors will appear? To answer this question, we note that there are 
2M hopping terms in Hkin, all in the form of ajaj. These hopping terms, when acting 
on a given basis vector, either annihilate it or change it into another basis vector with 
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some amplitude. In the latter case, the occupation numbers of the newly generated 
basis vector are readily obtained from those of \v). However, this information is not 
what we really need. The problem that really matters is, which basis vector is it among 
the basis vectors tabulated in the array A? Or more precisely, which row does it belong 
to in A? 

Here we will invoke the so-called hashing technique to fulfill this aim fl6\ \T7\ [18] . 
The basic idea is to define a tag for each basis vector, that is, to condense the information 
of the vector into a single entity. Thereafter, to see whether two vectors are the same, 
rather than comparing their elements one by one, we only need to see whether their 
tags are the same. Concretely, the tag of the f-th basis vector is defined by a function 
T, 

T{v) = T{A,i,A,2,...,A,m). (9) 

Numerically, this function should be readily evaluated. Moreover, since we want to 
identify the basis vectors with their tags, it is mandatory that different basis vectors 
have different tags. In other words, there should be a one-to-one mapping between the 
rows of A and the elements of the array T. 

A fortunate case is that the tag T(v) coincides with v. In this case, by calculating 
the tag of a basis vector, we know its rank among all the basis vectors. However, 
generally it is hard to find such a function. The compromise is to give up this hope 
and impose only the condition that all the tags are different, which is relatively easy to 
meet. A candidate of the tag function is 

M 

T{v) = J2 VPi^vi, (10) 
j=i 

with Pi being the i-th prime number. This function is linear in the occupation 
numbers and are readily calculated. More importantly, since the y^'s are radicals of 
distinct square- free numbers, they are linearly independent over the rationals [19], and 
therefore different vectors have different tags necessarily. An alternative tag function is 
^(^) = J2iLi(}^Pi)^m- By some simple number theory [20], it is ready to see that this 
tag function is also a viable one. The tag function the authors use is of the form flTU]) 
but with Pi = 100 * 2 + 3, which is easier to program. 

Given a basis vector \v) specified by a set of occupation numbers but with v 
unknown, we calculate its tag according to equation (fTII]) . then search the tag among 
the array T to locate its position, i.e. the value of v. Here another trick is possible. 
Originally, the array T is unsorted, i.e., the tags are not arranged in ascending or 
descending orders according to their values. To search a given element, the only way 
is to check the elements one by one and that will take on average D/2 trials to find 
out that element. That is a huge work since D can be on the order of 10®. A simple 
trick saves the workload significantly. Rather than searching inside an unsorted array, 
we had better search inside a sorted array so that we can make use of the Newton 
binary method [16]. That will take at most logg -D trials to find out the target. For 
D = 2^° ~ 1.05 X 10^ it takes at most 20 trials to locate the target [21]. 
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Figure 1. Sparsity pattern of the Hamiltonian matrix H in the case N = M ~ 10. 
Every spot corresponds to a non-zero element. The dimension of the Hilbert space is 
I? = 92 378 and the number of non-zero elements is nz = 1 064 777. On average, there 
are 11.5 non-zero elements per column. 



Thus we first sort T in ascending (or descending) order with the quicksort algorithm 
|22] . For clarity, we denote the sorted array as TSorted. In doing so, we can also 
prepare another D-element array ind which stores the positions of elements of TSorted 
in the original array T. More precisely, T(ind(i) )=TSorted(i) . In MATLAB, this can 
be done simply with the code 

[T,ind]=sort(T) 

Here we overwrite the original array T with TSorted. For those programming with 
Fortran, the ORDERPACK package by Olagnon can be used to fulfill the same aim 
[23]. 

We summarize the procedure to establish the matrix H_kin as follows. The non- 
zero elements are determined column by column. Given an arbitrary basis vector 
we apply the hopping terms a\aj onto it. If A^j > 1, we have 

a\aj\v) = sj [A^i + l)A^j\ . . . , + 1, . . . , - 1, . . .). (11) 

We then calculate the tag of the vector on the right hand side and search it among 
the sorted array T. Suppose T(w) = T^, we then know the resulting basis vector is the 
u = ind(w) -th one. We have thus found a non-zero element with coordinates {u, v) and 
value —J^ {^A^i -|- \)Ayy This process is repeated as v runs from 1 to D. Obviously, it 
can be parallelized. It is also clear that the overall time cost in this step is on the order 
of MD Xog^D. 

In figure [H we plot the sparsity pattern of the Hamiltonian H in the case N = M = 
10. We see that on average there are only 11.5 non-zero elements per column, which 
is four orders smaller than the dimension D = 92 378. We would like to mention that 
with our MATLAB code, it takes about 15 seconds to set up the Hamiltonian matrix 
and plot the pattern using the code spy(H). 
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5. Numerical results 

After preparing the Hamiltonian in a sparse matrix form, we can use the Lanczos 
algorithm to compute the ground state and low lying excited states and their energies. 
There are some well developed packages for this purpose and our philosophy is not to 
reinvent the wheel. For those programming with Fortran, the ARPACK package [21] by 
Lehoucq et al. is a very useful aid. For those programming with MATLAB, it is enough 
to invoke the "eigs" command. For instance, the code 

[Evec ,Eval] =eigs (H,2, ' sa' ) 

returns the two smallest eigenvalues (the ground state energy and the first excited state 
energy) of H in the 2x2 diagonal matrix Eval and their corresponding eigenvectors 
(the ground state and the first excited state) in the D x 2 matrix Evec. Here we would 
point out that when executing "eigs" , MATLAB invokes the very ARPACK package to 
do the job 

With the ground state on hand, we can then calculate various quantities to gain 
some physical insights of the model. One quantity that is of primary interest is the 
single-particle density matrix (SPDM) associated with the many-particle ground state. 
In the Wannier state basis, it is defined as 

p'lf = {G\ala,\G), (12) 

with 1 < j < M. All one-particle variables, e.g., the momentum distribution, are 
captured in the SPDM. 

In general, the SPDM is hermitian, semi-positive-definite, and of trace equal to the 
particle number. In the present case, the SPDM is subjected to more constraints. Due to 
the translation and reflection invariance of the ground state [I2] , we have pl^^ = pf^^ ^-^^ 
for an arbitrary k and also p^^ = p^^^ . Therefore, the SPDM is real, symmetric, and 
cyclic. These good properties reduce the number of matrix elements to be computed 
from M(M + l)/2 to [M/2] + 1, where [■] is the floor function. 

5.1. Condensate fraction 

According to the Penrose-Onsager criterion [26], a condensate is present if and only 
if the largest eigenvalue Ai of p*^^^ is macroscopic, i.e., fc = \i/N is on the order of 
unity and the ratio fc is called the condensate fraction. In the non-interacting case, 
all particles reside in the lowest Bloch state (a zero momentum state), and the system 
is in a pure condensate state with fc = 1. As the interaction is turned on, more and 
more particles will be kicked into higher Bloch states and the condensate is said to be 
depleted. In the thermodynamic limit (M goes to infinity with N/M = 1 fixed), there 
is a critical value [27] (~ 4.65) of U/ J beyond which fc vanishes. 

To gain a picture of this phase transition, we have numerically calculated the 
condensate fraction as a function of the ratio U/J, with three different lattice sizes. The 
results are shown in figure [2]^a). We see that as U/J increases, the condensate fraction 
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U/J 



Figure 2. (Colour online) (a) Condensate fraction fc and (b) correlation PQ||/y2] 
functions of the ratio U/J. In each panel, from up to down, the size of the system is 
N = M = 9, 11, and 13, respectively. 



decreases monotonically. However, the finite size effect is significant. The condensate 
fraction is far from being vanishing in the deep Mott insulator regime {U/J^ 4.65). 
Actually, since Ai > {Y,k ^k) /M = N/M, has a lower bound 1/M. 

5.2. Off-diagonal long range order 

The presence of a condensate is also associated with an off-diagonal long range order 
|28j . That is, a condensate is present if the off-diagonal element of the single particle 
density matrix Pj-j'* converges to a finite value as |« — j| — )■ oo. This is consistent with 
the Penrose-Onsager criterion. Actually, converting into the Bloch state representation, 
the density matrix takes the form 



P^L = (G\blb,,\G) 



M 

= (13) 

Here in the third line, we have used the cyclicity of p^^\ Thus the SPDM is diagonal in 
the Bloch state representation. Its eigenvalues coincide with its diagonal elements, and 
its eigenstates (called natural orbits) coincide with the Bloch states. It can be proven 
that all the elements of p*-^^ are non-negative [12]. Therefore, the largest eigenvalue of 
the SPDM is just p^q^q^ with qi = 0, and is of the explicit expression 

We then see immediately that, if p^^j decreases monotonically with j, fc and p^j 
converges to the same value in the thermodynamical limit. 

Thus the phase transition can also be investigated by examining the behavior of the 
off-diagonal elements. In figure [2]^b), we show how the element Po|}^^/2] behaves as U/J 
varies. We choose this element because it corresponds to the correlation between two 
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sites with the largest distance on a circle. Moreover, as M tends to infinity, the distance 
[M/2] also tends to infinity. Comparing with figure [2l^a), we see that the correlation 
Po[M/2] decreases much faster than fc, especially, in the deep Mott insulator regime, it 
does drop to zero. 

5.3. Occupation variance 

In the previous subsections, we see that both the condensate fraction and off-diagonal 
elements suffer from a strong finite size effect. Here, we show that the fluctuation of the 
occupation number on one site 

a, = ^J{G\h]\G)-{G\hi\Gy (15) 

is not very sensitive to the size of the lattice [29], as long as M >> 1. In flgure [3l we 
show function of U / J for flve different lattice sizes. The difference between the 

curves are hardly visible. This indicates that the curve has already converged to its 
value in the inflnite lattice limit. 

This suggests that the SF-MI transition can not be tracked in any local variables. 
This is why we do not fix the distance between the two sites when calculating the off- 
diagonal element in the proceeding subsection. By the Hellmann-Feynman theorem, we 
have dEc/dU = (dH/dU) = Maf/2. Our numerical result for aj suggests that the 
ground state energy is a smooth function of U. 

1.0 
0.8 
0.6 

b" 

0.4 
0.2 
0.0 

4 8 12 16 20 

U/J 

Figure 3. Variance of the occupation number ai at an arbitrary site i as a function 
of the ratio U/J [29]. Actually five different lattice sizes, i.e., N = M = 8,9, 10, 11, 12, 
are investigated. However, the curves all collapse onto the same one. 

6. Conclusion and discussion 

Exact diagonalization is simple conceptually but never trivial in programming. Taking 
the Bose-Hubbard model as a working example, we have illustrated the architecture of 
numerical exact diagonalization. Some essential tricks, namely, ordering, enumerating, 
hashing, sorting, and searching, were explained in detail. These tricks are believed to 
be quite general and can be adopted in many other situations |30] . 




Exact diagonalization: the Bose-Huhhard model as an example 



12 



For example, we show how the idea of ordering can save computation if we want 
to work within the subspace T/q- To impose the condition i^' = 0, we had better work 
with H in equation (jll) [T?]. This time, it is the interaction part 



H, 



int 



2M 



H^nt = ^Y.ll^\AKA. (16) 

Qim 93,94 

that costs most effort. Here the condition qi + q2 = i?3 + q'4 (mod M) is taken imphcitly. 
There are exactly terms in the sum. By noting hq^ = [6^^, = 0, we have 

Hint = Yl ^Ji92^9394, (17) 

9l>92 93>94 

where -B^j^j = (2 — '^9192)^91^92 defined for qi > q2. We can reduce the computation 
further by making use of the hermicity of Hmt- We define (^1^2) = Mqi + q2- The terms 
S^^q^ are ordered according to their tags {qiq2)- We then rewrite f lT7|) as 

X] ^Ji92^9394 + ^Ji92^9394 + ^Jig2^9394j • [^^) 

(9l92) = (9394) {gi92)>(9394) (9l92)<(9394) 

The third term is hermitian conjugate to the second one. Thus in setting up the matrix 
corresponding to Hmt, "we only need to consider the non-zero elements due to the second 
term, those due to the third term are then determined automatically. Overall, in the 
case of = M = 10, the number of terms need to be considered is reduced from 1000 
to 180. 

The readers are encouraged to convert the procedure described in this paper into 
codes and explore the interesting physics in the Bose- Hubbard model, which is surely far 
from being exhausted in the present paper. For example, we have shown how to study 
the condensate fraction as a function of the ratio U/ J. On this basis, an immediately 
accessible problem is then how the superfluidity density varies with U/J. The subtle 
relation between condensation and superfluidity can then be investigated. For more 
details, see [29]. We would like to mention that the coding does not cost much effort. It 
takes no more than 100 lines in MATLAB, and is very efficient. For the A^ = M = 12 
case, it takes around 4 minutes to set up the Hamiltonian matrix and 1 minute to solve 
the ground state on our desktop computer. Note that the dimension is D = 1352 078. 
By working in the subspace Ho, we have successfully performed exact diagonalization 
for a system as large as A^ = M = 13 on our computer. Systems with larger sizes may 
be investigated by working in the subspace. Therefore, the Bose-Hubbard model 
may serve as a good topic for teaching exact diagonalization in an undergraduate or 
graduate course of computational physics. 
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